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Abstract. We analyse the static solutions of attractive Bose-Einstein condensates 
under transverse confinement, both with and without axial confinement. By full 
numerical solution of the Gross-Pitaevskii equation and variational methods we map 
out the condensate solutions, their energetic properties, and their critical points for 
instability. With no axial confinement a bright solitary wave solution will tend to decay 
by dispersion unless the interaction energy is close to the critical value for collapse. In 
contrast, with axial confinement the only decay mechanism is collapse. The stability 
of a bright solitary wave solution increases with higher radial confinement. Finally we 
consider the stability of dynamical states containing up to four solitons and find good 
agreement with recent experiments. 
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The presence of attractive interactions in Bose-Einstein condensates (BECs) leads 
to rich and intriguing nonlinear phenomena. A key example is the formation of bright 
soliton-like structures P,[2l[3]. A bright soliton is a one-dimensional (ID) density wave 
that propagates without spreading due to a balance between attractive interactions and 
dispersion. In 3D and under transverse confinement, the analog is a bright solitary wave 
(BSW) solution, which is self-trapped in the axial direction [H [5], [6l [7]. Due to their 
self-trapped nature, BSWs hold significant advantages for atom-optical applications, 
such as atom interferometry. Another important property of attractive BECs in 3D is 
the collapse instability. In 3D a homogeneous condensate with attractive interactions 
is always unstable to collapse [8]. However, the presence of trapping can stabilise the 
condensate against collapse up to a critical number of atoms [9l|10]. Indeed, the collapse 
instability has been crucial in the experimental formation of BSWs P, E], [3] . A highly- 
populated repulsively-interacting BEC has its interactions switched to attractive via a 
Feshbach resonance. This induces the collapse of the condensate, with one [2] or more 
[U [3] BSWs emerging from the collapse. For the case of multiple BSWs, a 7r-phase 
difference leads to a repulsive solitonic interaction that is important in stabilising their 
collisions against collapse [71 [HI |12l [13] . 

Table 1 summarises the three experiments to date that have generated BSWs of 
attractive BECs, at Rice University [1], ENS in Paris [2], and JILA [3]. They feature 
cylindrically-symmetric traps with radial harmonic confinement of frequency Ur-, and 
either a confining or expulsive axial harmonic potential. Note that due to the presence 
of axial confinement these are not true solitonic states and from now on we will generally 
define BSWs to be solutions under zero axial confinement. The strength of the atomic 
interactions, characterised by the s-wave scattering length Os (os < for attractive 
interactions), relative to the trapping potential is crucial in determining the onset of 
collapse, and can be parameterised in terms of the dimensionless parameter [13], 

* = ^, (1) 

where = \Jh/ mojr is the radial harmonic oscillator length and m is the atomic 
mass. Collapse occurs when the atom number exceeds a critical population or, 
correspondingly, when the interaction parameter k exceeds a critical value kc- The 



Table 1. Parameters of the three matter wave 'sohtons' experiments to date. For the 
ENS experiment, we use an estimated atom number of TV = 4500, which faUs within 
the error of their measurements, while for the JILA experiment we assume N=1500. 
The trap ratio A is defined as A = uozl^r- 



Experiment 


Atomic species 


UJrl2-K (Hz) 


|A| 


Os (nm) 


N 


k 


Rice [I] 


^Li 


800 


0.005 (confining) 


-0.16 


5000 


0.6 


ENS [2] 


^Li 


710 


0.1 (expulsive) 


-0.21 


4500 


0.65 


JILA p] 


85Rb 


17.5 


0.4 (confining) 


-0.6 


1500 


0.35 
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value of kc is a crucial consideration and has been the subject of several theoretical 
investigations, both for condensates under a 3D potential [H Ull [15], [161 [13 [HI [IH] and 
for true BSWs under zero axial confinement P, [7j. It depends on the geometry of the 
trap and is of the order of unity. An accurate method to probe kc is to numerically solve 
the full 3D Gross-Pit aevskii equation that describes the BEG mean-field, and isolate 
the point where solutions no longer exist [6], [151 [13 [H]- This method is numerically- 
intensive but can be approximated using a variational approach [H |6]. Indeed, the 
variational method has been used to analyse solutions over the full range of axial trap 
geometries - confining, expulsive and zero axial trapping. However, this has not been 
directly compared with the more accurate method of solving the full Gross-Pitaevskii 
equation. In addition to these methods, the nonpolynomial Gross-Pitaevskii equation, 
an effective ID equation derived from the 3D Gross-Pitaevskii equation, has been used 
to predict the critical interaction parameter for a bright solitary wave [7] , while a simple 
analytic approach [19] has been used to derive the critical interaction parameter for an 
attractive condensate under 3D trapping. 

In this paper we theoretically analyse the solutions of an attractive BEG featuring 
radial confinement and under the full range of axial potentials - no axial trapping, 
confining axial potential and expulsive axial potential. This is performed by solving 
the full Gross-Pitaevskii equation and employing an approximate variational approach. 
We map the regimes of instability and the energetics of the solutions. In Section 1 we 
present our theoretical framework of the Gross-Pitaevskii equation and the variational 
approach. In Section 2 we use these methods to examine the stable BSW solutions in 
the absence of axial trapping and the occurrence of a finite excitation energy to the 
solution. In Section 3 we consider the solutions under an axial potential, which is either 
confining or expulsive, and compare to experiment. Finally in Section 4 we model the 
multiple interacting 'solitons' observed in the JILA experiment [3] and probe the critical 
number for such states, finding good agreement between experiment and theory. 



1. Theoretical framework 



1.1. Gross-Pitaevskii equation 



In the limit of zero temperature a dilute BEG can be approximated by a mean-field 
'wavefunction' ip{r,t) which satisfies the Gross-Pitaevskii equation |20j . 



dip 



.2 



h 

2m' '2 



ip. (2) 



Here is the number of atoms in the condensate and g = Anti^as/m parameterises 
the strength of the s-wave atomic interactions. The trap used to confine the BEG is 
assumed to be harmonic and cylindrically-symmetric, where uir is the transverse trap 
frequency and A = oj^/uj. is the trap ratio. The wavefunction has been normalised 
such that / |?/'(r)p(i'^r = 1. Furthermore the energy of the system is defined by the 



Bright solitary waves and trapped solutions of attractive BECs 



4 



Gross-Pitaevskii energy functional, 

E[i,] = I dh ||!-|V^(r)r + Imcol ^ + A^^^) |^(r)p + ^|^(r)r|(3) 

The basic Gross-Pitaevskii equation provides an excellent model of BECs at the mean- 
field level [21], and although it is insufficient to describe collapse dynamics, where 
higher order effects such as three-body loss become considerable, it provides a good 
model to infer the onset of collapse [151 IE]- For finite g there are no analytic solutions 
of the 3D Gross-Pitaevskii equation and we obtain the solutions numerically via the 
imaginary-time technique: by propagating the Gross-Pitaevskii equation in imaginary 
time {t —it), the system relaxes to the lowest energy state of the system [2T] . 



1.2. Variational method 

As well as solving the Gross-Pitaevskii equation numerically we use a variational 
technique to derive approximate solutions and give valuable insight without intensive 
numerics. Assuming tight radial trapping we decouple the 3D wavefunction ip{r) into 
the product of an axial and radial component, with the radial component assumed to 
be a gaussian harmonic oscillator state [H El El Ej. For the axial component there are 
two regimes: 

(i) For |A| ^ 1 the axial direction is dominated by interactions. In ID and the 
limiting case of A = 0, the Gross-Pitaevskii equation supports bright soliton solutions of 



the form ipi^) = \^sech.{z/Q [20]. The healing length ^ = h/ yniY)m\g\ characterises 
the soliton width, where uid is the peak ID density. Assuming our 3D solution to have 
this axial profile leads to the ansatz [3, El Ej , 



1/2 / ^2 



= 1,2^ ) ''''' [tJ [-^j- 

where and represent the axial and radial sizes, respectively. Substituting this ansatz 
into the energy functional of equation ([3]) leads to the ansatz energy [6] , 

=2^[ll^3l!)^ ['^ + 

(ii) For |A| > 1, the harmonic potential will dominate the axial direction. For a 
confining axial potential (A^ > 0), it is then more appropriate to assume a gaussian 
axial profile. This leads to, 

^.a.s(r, .) = exp exp , (6) 

and. 



Equations ([5]) and ([7]) define an energy landscape in terms of Iz and Ir- Minimising 
the energy with respect to these variational parameters leads to the variational solution. 



Bright solitary waves and trapped solutions of attractive BECs 5 




Figure 1. Energy landscape according to equation ([5]) for A = in the (a) stable 
regime k < {k = 0.35) and (b) unstable regime k > k^ {k = 0.9). Energy is in 
units of huir- White contours highlight the shape of the landscapes. In (a) two key 
regions are indicated: (I) the global energy minimum (collapse region) and (II) the 
local energy minimum corresponding to the variational solution. Whereas (a) features 
two minima, corresponding to collapse and the BSW solution, (b) features only a single 
minima representing collapse. In (a) the dashed green line indicates the low-energy 
path through the energy landscape from the origin to large Iz- For < this path 
represents the collapse channel (CC) and for Iz > 1^ it represents the dispersive channel 
(DC). 

2. Bright solitary wave solutions in the absence of axial trapping 

2.1. Stable and unstable energy landscapes 

We first examine BSW solutions (in the absence of axial trapping) using the variational 
ansatz of equation (j4]). This approach been performed previously by Carr and Castin 
[6], while Perez-Garcia et al. |1] performed a similar analysis using the gaussian ansatz 
of equation ([6]). A typical energy landscape according to equation ([5]) for a stable BSW 
solution is shown in figure [l](a), corresponding to k = 0.35. At the origin (region 1) the 
interaction term in equation ([5]) diverges to negative values since Og < for attractive 
condensates. This region represents the physical collapse of a wavepacket, although the 
global energy minimum itself is unphysical since it has zero width Importantly, there 
exists a stable local energy minimum indicated by region II in figure[T](a). This represents 
the BSW solution. A typical unstable energy landscape is shown in figure [D^b), for a 
large interaction parameter of /c = 0.9. No local energy minimum exists, and the whole 
parameter space is unstable to collapse. 

2.2. Critical point for collapse 

Using the sech-based ansatz of equation (j4]), Carr and Castin have shown that the critical 
interaction strength for collapse of a BSW is = 0.76 [6]. Similarly, Perez-Garcia et 

I In reality the collapse stops at a finite width when the three-body losses become considerable and 
reduce the atom number to below the critical number [21j . 
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Figure 2. (a) Evolution of the axial density (integrated over the radial direction) of 
a BSW as the interaction parameter k is adiabatically increased past the critical point 
for collapse fcc- The ramping rate is dk/dt = 0.12 s~^. (b) Axial and radial width of 
the BSW as a function of interaction strength k according to the 3D Gross-Pitaevskii 
equation (crosses) as it evolves under adiabatic ramping. The predictions from the 
sech-based variational method of equations (|4]) and (O are shown by solid lines. The 
critical interaction strength for collapse according to the sech-based variational method 
k^ and the Gross-Pitaevskii equation k^^ are indicated by vertical dotted and dashed 
lines, respectively. 



al. [1] have used the gaussian-based ansatz of equation to predict a critical point 
of fc^ = 0.778. Using the non-polynomial Gross-Pitaevskii equation, Salasnich et al. 
predict that = 2/3 [7]. These methods consider approximated solutions to the 3D 
Gross-Pitaevskii equation. A more exact method, although considerably more intensive, 
is to numerically solve the full 3D Gross-Pitaevskii equation [6| [T7]. 

In order to isolate kc from the full 3D Gross-Pitaevskii equation we employ an 
adiabatic ramping technique. Firstly, a stable BSW solution {k <^ kc) is obtained 
by propagating the 3D Gross-Pitaevskii equation in imaginary time. Then, in real 
time, the interaction parameter k is increased adiabatically. We employ a ramping rate 
dk/dt = 0.12 s~^. We have verified that a lower ramping rate does not change the 
point of collapse. An example of a BSW under adiabatic ramping is shown in figure 
|2]^a). For k < k^ the condensate progresses through the stationary state solutions, 
becoming progressively narrower and of higher density, but aX k = kc the condensate 
suddenly collapses. Using this technique we have isolated the critical interaction strength 
according to the full Gross-Pitaevskii equation for a BSW to be k^^ = (0.675 ± 0.005). 
The error margin arises from the finite time over which the wavepacket collapses. We 
have confirmed that if the value of k is ramped up to, and then maintained at, any 
value less than this k^^ , the condensate remains stable. This value is different to the 
value of k^^ = 0.627 obtained by Carr and Castin [6] through numerical relaxation of 
the 3D Gross-Pitaevskii equation. However, our result is consistent with Gammal et al. 
[T7] who employed numerical relaxation under a weak axial trap of A = 0.01 (very close 
to the axially-homogeneous case) and obtained k^^ = 0.676. Furthermore, this value is 
in good agreement with Salasnich et al. [7j who analytically derived kc = 2/3 from the 
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nonpolynomial Gross-Pitaevskii equation. 

In contrast to the adiabatic ramping method of quantifying /c^^, the numerical 
relaxation approach relies on searching for solutions to the Gross-Pitaevskii equation 
(by propagating in imaginary time) for increasing values of /c, and locating A;^^ as being 
the point when solutions can no longer be found. This has been employed successfully 
to obtain the critical interaction strength for attractive condensates under 3D trapping 
jl5[ [T7]. However, for BSWs close to the critical interaction strength, the local energy 
minimum representing the BSW solution can be so shallow (see solid line in figure [3]^c)) 
and localised that the imaginary time method is very sensitive to the initial guess and 
can become ineffective unless the initial guess is practically the solution itself. 

2.3. Lengthscales of the bright solitary wave 

The axial and radial widths of the BSW solutions, /° and /°, are shown in figure [2](b) 
as a function of k, from the full Gross-Pitaevskii equation (crosses) and the sech- 
based variational method (lines). We evaluate the axial width of the Gross-Pitaevskii 
equation solutions as being the lengthscale over which the axial density decreases to 
sech^(l)?7,o — 0.42no, where Uq is the central peak density, and the radial width as the 
lengthscale over which the radial density decreases to no/e. 

The variational and full Gross-Pitaevskii equation methods give the same 
qualitative description. The wave is elongated in the axial direction (/° > /,°), with 
radial width l^ remaining close to ar while the axial width /° diverges as k — > 0. 
For increasing interaction strength, 1^ approaches Ir, and at the point of collapse 
the BSW is approximately isotropic (/° = /°), as noted elsewhere [6]. However, 
the variational method consistently overestimates the widths predicted by the Gross- 
Pitaevskii equation. This implies that the variational solutions have lower peak density 
than the Gross-Pitaevskii equation solutions and explains why they predict a higher 
value of kc. 

2.4. BSW energetics 

Consider the variational BSW solution indicated by region II in figure [T]^a). We denote 
the energy minimum in this region by -Emin- In the radial direction, the local energy 
minimum is well-bounded due to the dominance of the radial confinement. In the axial 
direction, however, there is low energy path in the energy landscape which leads from the 
origin to large Iz, passing through the local BSW energy minimum. We have indicated 
this path on the energy landscape in figure [TJ^a) by the dashed green line. This low- 
energy path represents the most energetically-favourable states of the BSW at a given 
axial width. In figure [^a) and (b) we show the variation of the energy along this path, 
plotted on different scales to highlight contrasting features. For < /°, the energy 
(figure [3]( a)) initially increases with decreasing due to the growth of the kinetic terms 
in equation ([5]). However, for even smaller 1^ the negative interaction energy begins to 
dominate and causes the energy to decrease (and ultimately diverge to — 00 as — > 0). 
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Figure 3. Energetics of the BSW (for A = 0) according to the sech-based ansatz. (a) 
Energy profile along the low-energy path through the energy landscape (dashed line 
in figurc[TJa)). The interaction strength is fc = 0.35. The energy of the solution -Emin 
(dotted line), the energy barrier to the collapse region E'coi (dot-dashed line) and the 
resulting energy difference AE'coi are indicated. The position of the energy minimum 
is shown (grey line), (b) Same as (a) but plotted on a different scale to highlight the 
energy of the dispersive channel E'dis (dashed line) and the resulting energy difference 
AE'dis = Etiis — -Emin- (c) The variation of i?min (soHd hue), E'coi (dotted hne) and 
i?dis (dashed line) with the interaction parameter fc. (d) The excitation energy AE of 
the solution, i.e. the lower of either the collapse excitation energy AE'coi or dispersive 
excitation energy AE'dis- 



Importantly, this gives rise to a finite energy barrier which stabihses the system against 
collapse. For Iz > l^, the energy (figure [3](b)) increases weakly with 1^. This is because 
interaction energy in equation ([5]) falls off more slowly than the kinetic energy. This 
energy profile means that if one creates a state with > then the interactions are 
sufficient to overcome dispersion and the state will oscillate in the energy minimum, 
with its length oscillating around 

In the limit l^ oo, the axial kinetic and interaction energy terms in equation ([5]) 
tend towards zero and the energy landscape becomes EscchilrJz co) = h? /{2mll) -|- 
muj'^l'l/2, with a minimum energy E^is = at = a^. As a result we see that 
the low energy channel in figure [3](b) tends asymptotically towards E^^g as oo. 
The fact that the energy is finite at = oo has important consequences. If one adds 
additional axial kinetic energy to the system, e.g. by exciting an axial breathing mode. 
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we can form states with higher energy than described by equation ([S]), i.e. states that he 
above the hne in figure [3](a) and (b). Importantly, if the energy of this excited state is 
greater than E^is then it has sufficient energy to overcome the interactions and expand 
indefinitely. This means that there is a finite excitation energy AE'^is = -Edis — -E'min to 
induce dispersion of the BSW, as indicated in figure [^b). We term this the dispersive 
channel (DC), as illustrated in figure [Tt^a). 

There is another channel by which the BSW can be excited to a non-solitonic state. 
For axial widths less than there is a high-energy saddle point in the energy surface 
that links the local BSW energy minimum to the collapse region. We term this the 
collapse channel (CC) and have indicated it in figure [Dj^a). In figure [3t^a) we see this 
feature as an energy barrier. The energy maximum has amplitude -Ecoi and so presents 
an energy barrier of magnitude Ai?coi = Eco\ — -E'min against the collapse of the BSW. 

In figure [3](c) we show the variation of these energies as a function of the interaction 
parameter k. The energy minimum -Emin decreases slightly with increasing k (solid line). 
The energy of the dispersive channel is independent of k (dashed line). The collapse 
channel energy Ef.o\ (dotted line) decreases rapidly with increasing k. When the critical 
interaction strength kc is reached, AE'coi = (there is no longer an energy barrier to 
prevent collapse), and for k > k^ no variational BSW solutions exist. 

These channels introduce a finite excitation energy to excite the BSW to a non- 
solitary state. The minimal excitation energy (for either of these channels), AE, 
characterises the stability of the BSW. This excitation energy is plotted in figure [3](d) 
as a function of k. For /c<0.66, the dispersive channel is of lowest energy, with 
AE increasing with k. If the BSW were excited in this regime it would be prone 
to dispersion. However, for A; > 0.66 the collapse channel is of lowest energy and AE 
rapidly decreases to zero as k ^ k^. A BSW excited in this regime will be prone 
to collapse. The maximum excitation energy AE ^ O.OQ/icJr is finite and occurs at 
k ~ 0.66. Consequently, to maximize the thermal and dynamical stability of a BSW 
requires tight radial confinement and k ~ 0.66. 

According to the JILA parameters (but with A = 0), the weak radial trapping 
leads to a small excitation energy AE ^ 0.08 nK, whereas the strong radial trapping 
of the Rice and ENS parameters leads to a much larger value of AE ^ 3 nK. In the 
experiments, following the collapse the resulting solitonic states are sufficiently dilute 
that it is unlikely that the system reaches thermal equilibrium within the experimental 
time frame. However, we can estimate the thermal energy of the system based on the 
temperature prior to collapse. In the JILA experiment there are initially approximately 
15000 atoms in the BEC and 500 thermal atoms, corresponding to a temperature of 
4.6 nK. This suggests that a BSW in the JILA system will be unstable to thermal 
effects, whereas the Rice and ENS systems may just support thermodynamically-stable 
BSWs. However, BSWs may have enhanced thermodynamic stability due to self-cooling 
via the radiation of hot atoms in the untrapped direction [12] . Of course, the experiments 
themselves featured an axial potential, unlike the A = limit considered here, and we 
will see in the next section that this modifies the thermodynamic stability dramatically. 
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Figure 4. Energy landscapes of equation (O for attractive BECs under the confining 
and expulsive axial potentials used experimentally: (a) JILA parameters with k ~ 0.35 
and confining axial trapping A = 0.4; (b) Rice parameters with k = 0.6 and confining 
axial trapping A = 0.005. (c) ENS geometry with k = 0.65 and an expulsive axial 
potential (A^ < 0) with trap ratio |A| = O.f. Energy is presented in units of hujr and 
white contours highlight the shape of the landscapes. The dashed line in each plot 
indicates the low energy pathway through the energy landscape, from the origin to 
large l^. In (a) and (b) the preferred decay route is the collapse channel (CC) while in 
(c) it is the expansive channel (EC). 



3. Solutions under an inhomogeneous axial potential 

3.1. Confining axial potential 

We now consider the condensate solutions in the presence of a confining axial potential 
(A^ > 0). In figure |4l^a)-(b) we plot the energy landscapes according to equation ([5]) for 
the parameters employed in the JILA and Rice experiments. In each case, the lowest 
energy path through the landscape is illustrated by dashed lines. In figure [5] we plot 
the energy profile along the low-energy path through the energy landscape. Due to the 
axial confinement the A^/^-term in equation ([5]) dominates for large 1^ and ultimately 
leads to the energy increasing with 1^. Consequently there is no dispersive channel in this 
system, as evident from figure [5] (dashed line). The presence of even small axial trapping 
therefore has a profound increase in the stability of the state. The large trap ratio of the 
JILA experiment (A = 0.4) significantly modifies the energy landscape from the A = 
case (figure [U^a)), whereas the weak trap ratio of the Rice experiment (A = 0.005) only 
has a weak effect on the energy landscape (on the lengthscales shown). 

It is still possible to excite the solution out of the energy basin into the collapse 
region, as indicated by arrows in figure |4](a) and (b). According to the sech-based 
variational method, this excitation energy is AE = 1.32huJr ~ 1 nK for the JILA 
experimental parameters and AE = 0.2AhuJr ~ 9 nK for the Rice experimental 
parameters, which are of the order of the typical thermal energy of the BEC. These 
values are significantly larger than the A = case and indicate enhanced stability under 
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Figure 5. (a) Energy profile along the low-energy patii tiirough the energy landscape 
(dashed lines in figure |4l^a)-(c)) for the JILA system (dashed line), Rice system (solid 
line) and ENS system (dotted line), (b) Same as (a) but on a different scale to highlight 
the small-scale features for the Rice and ENS system. Although not evident on the 
lengthscales shown, the energy in the Rice system tends to infinity as — > oo due 
to the effect of the confining potential. Energy is presented in units of htOr and each 
experiment features a different ujr- Although the collapse barrier in the Rice system 
appears smaller than for the JILA system, it is in reality much larger due to the much 
larger value of uor employed in the Rice experiment. 



an axial trap. 

For an attractive BEC in a spherically-symmetric trap the critical interaction 
strength has been predicted to be kc ~ 0.57 by a variety of methods [16], [171 ITS] . 
Using imaginary-time relaxation of the Gross-Pitaevskii equation, Gammal et al. [17] 
predicted A;^^ for various trap ratios and, for the JILA trap geometry, give excellent 
agreement with the experimentally measured value of kc [TOl 122] . By simulating the 
Gross-Pitaevskii equation under the adiabatic ramping technique we have obtained the 
critical interaction strength for collapse fee, with the results presented in figure EJ^a) 
(crosses). Note that the shaded region represents stable solutions. For a confining axial 
potential our Gross-Pitaevskii equation results are consistent with the results of Gammal 
et al. [17] (circles). We also present the results of the sech-based variational method 
(dotted line) and the gaussian-based variational method (dashed line). Although these 
predictions give the same qualitative features, namely that kc decreases weakly with 
increasing trap ratio A, these variational methods consistently overestimate the critical 
interaction strength by around 20%. 

In figure ^ih) we plot the axial width of the solutions as a function of the trap 
ratio A for a fixed interaction strength oik = 0.35. The Gross-Pitaevskii equation (solid 
lines) and variational predictions (dotted and dashed lines) show good quantitative 
agreement. The radial widths /°, shown in the inset, also show reasonable agreement. 
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Figure 6. (a) Critical interaction strengths as a function of |A| for a confining potential 
> (red data, right-hand side of figure) and expulsive potential < (blue 
data, left-hand side of figure) . The region of stable solutions, according to the Gross- 
Pitaevskii equation, is shaded. For the confining potential there is only single critical 
point, namely fee for collapse, while for the expulsive trap there is additional critical 
point for expansion k^. Our Gross-Pitacvskii equation results are shown by crosses, the 
sech-based variational method by the dotted line and the gaussian-based variational 
method by the dashed line. The predictions of Gammal et al. [17] for a confining axial 
trap are shown (circles), (b) Axial size of the condensate solutions as a function of 
|A| according to the Gross-Pitaevskii equation (solid line), the sech-based variational 
method (dotted line) and the gaussian-based variational method (dashed line). We 
employ an interaction strength oi k = 0.35 such that wc probe the range of solutions 
shown by the horizontal solid line in (a). The insets show how the radial width depends 
on |A|. Note that the scale for |A| on the left and right sides of (a) and (b) are different. 



3.2. Expulsive axial potential 

The ENS experiment employed an expulsive (A^ < 0) trap in the axial direction [2]. 
This geometry leads to an additional effect whereby the wavepacket can be ripped apart 
by the expulsive trap [21 [6]. This occurs when k is less than a /ower critical interaction 
strength for expansion fcg. The upper and lower critical interaction strengths for an 
expulsive trap have been mapped out using the sech-based variational method [2], [6]. 
The energy landscape for the ENS system is shown in figure |l](c). The energy profile 
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along the low-energy path through the energy landscape is shown in figure (dotted 
line). The expulsive axial potential leads to a downward slope in the energy landscape for 
large Iz- In this region wavepackets will expand indefinitely. We term this the expansive 
channel, as indicated in figure 111(c). Note that we define the effect of expansion to be 
induced by the axial potential and so it is distinct from the effect of dispersion discussed 
in Section [21 The energy barrier between the solution and the expansive region decreases 
with decreasing interaction strength, such that below /Cg no solutions are stable against 
expansion. Note that, according to the energy landscape (figure Hl^c)), the minimum 
excitation energy is AE ^ O.OOG^cUr ~ 0.2 nK and corresponds to excitation into the 
expansive region. Although this thermodynamic energy is less than the typical energy 
of a BEC, the expulsive trap is expected to lead to the removal of hot atoms and so 
may lead to anomalously low condensate temperatures. 

We have calculated the critical interaction strengths for collapse kc and expansion 
ke using the Gross-Pitaevskii equation. In the former case we adiabatically ramp up k at 
fixed |A| until collapse in induced, and in the latter case we adiabatically ramp up |A| at 
fixed k until expansion is observed. In figure Et^a) (left-hand side) we plot kf. and k^ as a 
function of trap ratio |A| using the full Gross-Pitaevskii equation (crosses), and the sech- 
based (dotted line) and gaussian-based (dashed line) variational methods. Note that 
the range of trap ratios we present is an order of magnitude smaller than for the case 
of the confining axial potential (right-hand side of figure Et^a)). The qualitative features 
are that k^ (upper data) increases weakly with trap ratio |A|, while k^ (lower data) 
increases sharply from zero at |A| =0. The region between /cc(|A|) and fce(|A|) (shaded 
region in figure [6]^a)) represents stable solutions, the region above A;c(|A|) represents 
interaction- induced collapse, and the region below A;e(|A|) represents potential- induced 
expansion. At some critical trap ratio |A|c, k^ = kg, and for |A| > |A|c there are no 
stable solutions. According to the variational methods, |A|^ ~ 0.23 while according to 
the Gross-Pitaevskii equation |A|^^ ~ 0.15. 

4. Dynamical multi-soliton states 

In the JILA experiment [3] up to six long-lived localised 'solitons' [§| were observed. 
These wavepackets oscillated axially in a robust manner, repeatedly colliding with each 
other. Interestingly, the total number of atoms could greatly exceed the critical number 
for collapse and yet the dynamics remained stable. This is because each wavepacket 
contained less than A'c atoms and so was individually stable to collapse, and furthermore 
a TT-phase difference existed between adjacent wavepackets to prevent them coalescing 
and collapsing [6], El [HI [13l [23] . A 7r-phase difference has also been inferred to exist in 
the train of solitons observed in the Rice experiment P,[71[23]. These 7r-phase differences 
are thought to arise through the imprints of quantum fluctuations [23] or the emergence 

§ The experimental system is 3D and features axial confinement, and so the wavepackets are not solitons 
in the strict sense. However, in this section we will follow the experimental protocol in referring to 
these wavepackets loosely as solitons. 
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of this stable configuration through repeated instabihties [HI [23] 

The JILA experiment examined configurations of up to four distinct sohton 
structures. We will probe the stability of configurations of several coexisting solitons 
in the specific JILA experimental set-up, namely a ®^Rb BEC in a three-dimensional 
trap geometry tu^ = 27r x 6.8 Hz and cu^ = 27r x 17.5 Hz. Crucially we will assume a 
TT-phase difference between adjacent solitons. A more rigorous approach to this problem 
would be to study excited nonlinear states of the system rather than multiple solitons, 
as has been performed in quasi-lD limit using Hermite functions by Michinel et al. |25j . 
Multi-soliton states with alternating phase have also been shown to be supported by 
the nonlinear Schrodinger equation elsewhere [26] . 

4-1. Generation of dynamical multi-soliton states 

We impose vr-phase differences to separate the attractive BEC into several distinct 
wavepackets, and for each configuration examine the critical number of atoms. We 
define the critical number for q distinct wavepackets to be N^c^. We generate the g = 1 
solutions by imaginary time propagation of the Gross-Pitaevskii equation with constant 
global phase and isolate the critical number N^i as when solutions can no longer be 
found. In a similar manner, we form a g = 2 state by solving in imaginary time subject 
to a TT-phase difference at 2 = 0, and infer the critical number A'^c2- 

Generation of g > 2 states is more complicated. To form a g = 3 solution we employ 
two TT-phase steps aX z = ±L/2, as illustrated in figure [71|^a), to form three distinct 
regions z < — L/2, —L/2 < z < +L/2 and z > +L/2. Imaginary-time propagation of 
the Gross-Pitaevskii equation under these phase constraints can preferentially populate 
the central region to point of collapse. Instead we generate the q = 3 states dynamically, 
with an example shown in figure [7](c). Initially we form the ground state with constant 
phase and repulsive interactions. Then, in real time, the interactions are switched 
to attractive while the stepped phase distribution (figure [Zt^a)) is imprinted onto the 
condensate for a finite time (a time of At = 80ms works effectively). This forms three 
wavepackets with no atom transfer between them (figure [TJ^c)). This dynamical method 
of generating solitons is loosely analogous to the experimental generation of solitons. 
However, whereas the alternating phases arise naturally in the experiment, we must 
introduce them numerically. 

When the phase imprinting is terminated, and providing the critical number Ncs 
is not exceeded, the solitons oscillate axially and repeatedly collide (figure [7](c)), 
as observed experimentally. The alternating phase structure prevents overlap and 
suppresses the collapse instability. However, the solitonic interaction excite collective 
modes in the solitons. When N exceeds a critical value, the system is unstable to 
collapse and the solitons are destroyed. This critical value depends on the lengthscale 
L, and so to obtain N^s we must vary L to obtain the most stable configuration. We 
find that the most stable configurations arise when L ^ 2,^. In order to generate a 
g = 4 state we perform a similar process by imposing three TT-phase steps at z = and 
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Figure 7. (a)-(b) The stepped phase profiles, with characteristic size L, used to 
dynamically generate g = 3 and 9 = 4 states. The phase profiles are enforced until 
t = 80ms. (c) Evolution of a g = 3 and (d) g = 4 state with N = 3000, and L = 3.5/im. 
Following the JILA experiment, we employ as = —0.6 nm. The plots show the evolution 
of the axial density, integrated over the radial direction, (e) Stability diagram oi N/Nd 
showing stable configurations of up g = 1, 2, 3 and 4 states. The line (N/Nc) is plotted 
for comparison (dashed line). 



z = ±L, as illustrated in figure [7|(b). Under stable parameters, four localised solitons 
are formed, which subsequently oscillated in the system in a stable manner. An example 
is shown in figure [7](d). 

4-2. Stability of dynamical multi-soliton states 

The stability of states with g < 4 and a fixed scattering length = — 0.6nm is shown in 
figure [Tl^e). If each soliton was independent with critical number A'ci, the total critical 
number for q solitons would be qNd. In figure [Tl^e) we have plotted this estimate as 
the ratio N/N^ (dashed line). However, the numerical data gradually deviates from this 
line and N^q < qNd, indicating that the solitons are increasingly affected by each other 
as their number increases. A similar trend was predicted in the quasi-lD approach of 
Michinel et al. [25]. 

An important question is how many solitons are favoured to exist for a given number 
of atoms A^. The addition of a soliton involves adding a vr-phase slip to the system, which 
is energetically costly. We can therefore expect that the most energetically-favourable 
number of solitons will be the minimum number of solitons/phase slips that can support 
A^ atoms without being unstable to collapse. This data is plotted in figure M^a) for the 
three scattering lengths employed in the JILA experiment. All three data sets diverge 
away from the ratio N/N^. as q increases. The divergence is greatest for low scattering 
lengths, which is related to the fact that these wavepackets are widest (since the soliton 
size is characterised by ,^ oc and therefore interact most with their neighbours. 

The region above the points represents configurations of the system which are completely 
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Figure 8. (a) The ratio N/Nd versus tlie most stable number of solitons q (up to 
5 = 4) for scattering lengths ag = — 0.6nm (red pluses/dashed line), — 0.77nm (blue 
crosses/dotted line), and — 1.59nm (green circles/solid line), (b) The experimental 
remnant data of (blue points with error bars) and the corresponding numerical 
data where N — iVoxp(black circles). Numerical data corresponding to = 0.9iVcxp is 
also shown (red crosses). Note that the points are for different scattering lengths. 



unstable to collapse. The region below the lines in figure E^a) also represents stable 
configurations of the system, but these configurations consist of a higher number of 
solitons/phase slips, and therefore correspond to higher energy. 

4-3. Comparison with the JILA experiment 

The JILA experiment [3] measured the observed number of solitons q and the total 
number of atoms N in each case, and the experimental data is presented in figure [H](b) 
(points with error bars). Note that both q and were averaged over many observations, 
and generally feature an error representing the spread in their observations. Although 
this data shows the same trend as our numerical data in figure El^a) (for fixed Os), it 
should be note here that each point corresponds to a different scattering length. We 
have therefore plotted the corresponding numerical results in figure [H](b). However, it 
is important to note the following: the number of atoms measured experimentally in 
the system A'exp will typically include non-condensed atoms, whereas in our analysis A^ 
describes the total number of atoms in the BEC. In other words, these experimental 
and theoretical parameters are not necessarily equal. 

Initially we assume that all the experimentally observed atoms are in the BEC, 
i.e. A^ = Aexp (black circles in figure [8](b)). The numerical data is similar to 
the experimental data with most numerical points lying within the spread of the 
experimental observations. However for two data points the numerical data over- 
estimates g, e.g. where the experiment observes close to one and two solitons we 
predict two and three solitons, respectively. However, assuming that 10% of the 
atoms observed experimentally represent non-condensed atoms, i.e. A^ = O.QA'exp (red 
crosses in figure [Hl^b)), we find that these two anomalous points are corrected and every 
numerical point is consistent with the experimental observations. This is a conservative 
estimate for the number of non-condensed atoms given the 'hot' state of the condensate 
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in the experiment. 

Note that if there was zero phase difference between the sohtons they would overlap 
upon collision, and the critical number would then be A^ci for q>l. 

5. Conclusions 

In this paper we have considered the solutions of an attractively-interacting Bose- 
Einstein condensate in the mean-field limit using the full Gross-Pitaevskii equation 
and variational techniques. Radial confinement is assumed throughout, and in the axial 
direction we consider the absence of trapping and the presence of confining or expulsive 
harmonic potentials. In particular we evaluate the critical interaction parameter kc for 
condensate collapse across the whole range of axial geometries using the full Gross- 
Pitaevskii equation. 

We note that the variational method is not particularly sensitive to whether the 
sech-based or gaussian-based ansatz is used, even in the limits of zero and strong axial 
trapping. We find that the variational method gives good qualitative agreement with 
the full Gross-Pitaevskii equation, but that significant quantitative differences can exist, 
particularly in the important critical points for collapse and expansion. 

In the absence of axial trapping, bright solitary wave solutions exist, which are 
self-trapped in the axial direction by the attractive interactions. According to the full 
solution of the Gross-Pitaevskii equation, we find that bright solitary wave solutions are 
stable to collapse only when the interaction parameter k < 0.675ih0.005. Analysis of the 
energy 'landscapes' reveals a finite excitation energy to excite the bright solitary wave 
to non-solitonic states, characterised by either collapse or dispersion of the wavepacket. 
Unless the interaction parameter is close to the collapse condition, i.e. k > 0.66, the 
BSW will tend to decay by dispersion. The excitation energy of a bright solitary wave is 
proportional to the radial trap frequency, and for sufficiently weak radial trapping, the 
excitation energy can be less than the typical thermal energy, indicating thermodynamic 
instability of such states. 

In the presence of either a confining or expulsive axial potential, kc becomes 
dependent on the trap geometry. For the confining potential, we regain the results 
of Gammal et al. [17] and good agreement with the experimental measurement of kc 
p!Ol [22] . Under an expulsive axial trap, there is also a /ower critical interaction strength 
fee below which the trap leads to expansion of the wavepacket. Therefore kc and kg define 
the upper and lower boundaries of a range of solutions. However, when the trap ratio 
|A| exceeds a critical value |A|c ~ 0.15 (based on the full Gross-Pitaevskii equation), 
no solutions exist for any interaction strength. Whereas for expulsive potentials and 
no confinement the main decay channel is dispersion, with axial confinement the only 
decay mechanism is collapse. Note that, for a weak expulsive potential, one can tune 
the system such that the energy barrier for collapse and dispersion are equal. 

Finally, based on the JILA experiment [3] we have generated dynamical states of 
multiple 'solitons' with alternating phase. We have mapped out the critical numbers 
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for these states and find good agreement with the quantitative experimental data, 
confirming the important role of vr-phase differences in systems of multiple solitons. 
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